Plots for Figures 3 and 4
Rahul Subramanian
7/8/2020
This document contains the code used to generate the panels in Figures 3 and 4.
The “big b_a” and “small b_a” parameters identified here were used to generate trajectories for the validation analysis.
The combined figure for Figure 3 was plotted in Matlab using the data generated in “R_0_Surface_Plot_data_for_revision.csv”. The Matlab code is included in the file “Surface_Plot_Panel_Script_add_colorbar_2D.m” in the repo.
library(ggplot2)
library(plotly)## Warning: package 'plotly' was built under R version 3.5.2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutsource("load_libraries_essential.R")## Warning: package 'dplyr' was built under R version 3.5.2library(pomp)## Warning: package 'pomp' was built under R version 3.5.2## Welcome to pomp version 2!
## For information on upgrading your pomp version < 2 code, see the
## 'pomp version 2 upgrade guide' at https://kingaa.github.io/pomp/.source("rahul_theme.R")Surface plot of b_a vs b_p vs R_0, color by b_a
model_name = "N_12"
load(file = paste0("../Generated_Data/Profiles/", model_name, "_Model/top_2_LL_data/b_a_profile_antibody_surface_plot_data.RData"))
load(file = paste0("../Generated_Data/Profiles/", model_name, "_Model/top_2_LL_data/b_a_profile_antibody_surface_plot_data_with_outlier.RData"))
#head(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_outlier)
antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers =
antibody_top_2_LL_from_b_a_profile_top_2_LL_no_outlier %>%
filter(b_p > 0.02)
f1 <- list(
family = "Arial, sans-serif",
size = 18,
color = "black"
)
f2 <- list(
family = "Old Standard TT, serif",
size = 15,
color = "black"
)
R_0_surf_plot_data = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers %>%dplyr::select(R_0, b_a, b_p)
write.csv(R_0_surf_plot_data,
"../Generated_Data/Profiles/N_12_Model/top_2_LL_data/R_0_Surface_Plot_Data_for_Revision.csv",
row.names = FALSE)
fig <- plot_ly(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
x = ~b_a, y = ~b_p, z = ~R_0, color = ~R_0)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = ' b_a',
zeroline = TRUE,
range = c(0,1),
zerolinecolor = toRGB("black"),
titlefont = f1,
tickfont = f2,
zerolinewidth = 6),
yaxis = list(title = ' b_p',
range = c(0,1),
zeroline = TRUE,
showline = TRUE,
zerolinecolor = toRGB("black"),
titlefont = f1,
tickfont = f2,
zerolinewidth = 6),
zaxis = list(title = 'R_0 ',
range = c(0,8.2),
zeroline = TRUE,
zerolinecolor = toRGB("black"),
titlefont = f1,
tickfont = f2,
zerolinewidth = 6)))
figfig <- fig %>% layout(scene = list(xaxis = list(title = ' b_a',
zeroline = TRUE,
zerolinecolor = toRGB("black"),
titlefont = f1,
tickfont = f2,
zerolinewidth = 6),
yaxis = list(title = ' b_p',
zeroline = TRUE,
showline = TRUE,
zerolinecolor = toRGB("black"),
titlefont = f1,
tickfont = f2,
zerolinewidth = 6),
zaxis = list(title = 'R_0 ',
zeroline = TRUE,
zerolinecolor = toRGB("black"),
titlefont = f1,
tickfont = f2,
zerolinewidth = 6)))
figSurface plot of b_a vs b_p vs R_0_NGM, color by b_a
R_0_NGM_surf_plot_data = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers %>%dplyr::select(R_0_NGM, b_a, b_p)
write.csv(R_0_NGM_surf_plot_data,
"../Generated_Data/Profiles/N_12_Model/top_2_LL_data/R_0_NGM_Surface_Plot_Data_for_Revision.csv",
row.names = FALSE)
fig <- plot_ly(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
x = ~b_a, y = ~b_p, z = ~R_0_NGM, color = ~R_0_NGM)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = ' b_a',
zeroline = TRUE,
range = c(0,1),
zerolinecolor = toRGB("black"),
titlefont = f1,
tickfont = f2,
zerolinewidth = 6),
yaxis = list(title = ' b_p',
range = c(0,1),
zeroline = TRUE,
showline = TRUE,
zerolinecolor = toRGB("black"),
titlefont = f1,
tickfont = f2,
zerolinewidth = 6),
zaxis = list(title = 'R_0_NGM ',
range = c(0,5),
zeroline = TRUE,
zerolinecolor = toRGB("black"),
titlefont = f1,
tickfont = f2,
zerolinewidth = 6)))
figfig <- plot_ly(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
x = ~b_a, y = ~b_p, z = ~R_0_NGM, color = ~R_0_NGM)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = ' b_a'),
yaxis = list(title = ' b_p'),
zaxis = list(title = 'R_0_NGM ')))
figSurface plot of R_0 vs R_0_NGM vs b_p, color by b_a
fig <- plot_ly(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
x = ~R_0, y = ~R_0_NGM, z = ~b_p, color = ~b_a)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = ' R_0'),
yaxis = list(title = ' R_0_NGM'),
zaxis = list(title = 'b_p ')))
figfig <- plot_ly(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
x = ~b_a, y = ~b_p, z = ~R_0, color = ~R_0_NGM)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = ' b_a',
zeroline = TRUE,
range = c(0,1),
zerolinecolor = toRGB("black"),
titlefont = f1,
tickfont = f2,
zerolinewidth = 6),
yaxis = list(title = ' b_p',
range = c(0,1),
zeroline = TRUE,
showline = TRUE,
zerolinecolor = toRGB("black"),
titlefont = f1,
tickfont = f2,
zerolinewidth = 6),
zaxis = list(title = 'R_0_NGM ',
range = c(0,5),
zeroline = TRUE,
zerolinecolor = toRGB("black"),
titlefont = f1,
tickfont = f2,
zerolinewidth = 6)))p = ggplot(data = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
aes(x = R_0,
y = R_0_NGM,
color = b_a)) + rahul_man_figure_theme +
scale_color_viridis_c(name = "Strength \n of Asymp\n Transmission") +
geom_point() +
xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number")
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a.png")
print(p)
dev.off()## quartz_off_screen
## 2Isolate MLE for each scenario
head(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers)## b_a M_0 V_0 K_0 R_0 b_q b_p p_S p_H_cond_S
## 1 0.00000000 5 13 14 7.061999 0.2401933 0.7957984 0.1574275 0.1749811
## 2 0.00000000 5 13 14 6.215237 0.2395227 0.9233632 0.1740948 0.1532892
## 3 0.06896552 5 13 14 5.984521 0.2204782 0.8913949 0.1735879 0.1757053
## 4 0.06896552 5 13 14 6.098822 0.2343746 0.9402490 0.1488610 0.1569547
## 5 0.10344828 5 13 14 6.845525 0.2237219 0.7932278 0.1296082 0.1818042
## 6 0.13793103 5 13 14 5.905499 0.2168981 0.8785885 0.1494229 0.1631618
## phi_E phi_U phi_S h_V gamma N_0 E_0 z_0 C_0
## 1 1.09 1.09 0.2 0.125 58.078038 8e+06 65553.63 10267.580 0
## 2 1.09 1.09 0.2 0.125 8.978508 8e+06 59637.11 11293.462 0
## 3 1.09 1.09 0.2 0.125 35.301204 8e+06 56854.62 9411.456 0
## 4 1.09 1.09 0.2 0.125 6.333218 8e+06 63566.34 13443.034 0
## 5 1.09 1.09 0.2 0.125 48.705248 8e+06 71702.50 12381.941 0
## 6 1.09 1.09 0.2 0.125 17.720241 8e+06 63878.93 12506.505 0
## social_distancing_start_time quarantine_start_time PCR_sens sigma_M
## 1 17 22 0.9 0.2764235
## 2 17 22 0.9 0.2750787
## 3 17 22 0.9 0.2768054
## 4 17 22 0.9 0.2709327
## 5 17 22 0.9 0.2752184
## 6 17 22 0.9 0.2751406
## beta_w_3 beta_w_2 beta_w_1 beta_w_0 g_0 g_F sigma_epsilon
## 1 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 2 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 3 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 4 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 5 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 6 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## G_w_y_scaling msg iter_num param_index loglik nfail trace_num
## 1 0.162 mif1 1 10 -627.3202 NA NA
## 2 0.162 mif1 1 11 -627.2093 NA NA
## 3 0.162 mif1 2 43 -628.1031 NA NA
## 4 0.162 mif1 2 44 -627.4300 NA NA
## 5 0.162 mif1 1 60 -628.1314 NA NA
## 6 0.162 mif1 2 76 -628.5305 NA NA
## loglist.se Antibody_Mean_LL Antibody_LL_SE Median_Herd_Immunity
## 1 0.008206760 -24.75616 0.004776843 0.1933564
## 2 0.010463251 -26.12591 0.007764530 0.1749998
## 3 0.007808338 -26.08064 0.008118650 0.1775874
## 4 0.008839892 -24.41689 0.002390176 0.1998329
## 5 0.010838884 -25.56156 0.006393602 0.2380878
## 6 0.008314098 -24.34239 0.002107043 0.2090317
## combo_num sim_subset_index duration_of_symp_1 duration_of_symp_2
## 1 3 1 5 0.01721821
## 2 4 1 5 0.11137708
## 3 11 1 5 0.02832765
## 4 13 1 5 0.15789761
## 5 15 1 5 0.02053167
## 6 19 1 5 0.05643264
## duration_of_symp gamma_total Beta R_0_P R_0_A R_0_S_1
## 1 5.017218 0.1993136 1.407553 1.0276405 0.0000000 1.1079374
## 2 5.111377 0.1956420 1.215961 1.0300678 0.0000000 1.0584630
## 3 5.028328 0.1988733 1.190161 0.9733061 0.3391599 1.0329881
## 4 5.157898 0.1938774 1.182424 1.0199752 0.3470370 0.8800840
## 5 5.020532 0.1991821 1.363506 0.9922669 0.6138540 0.8836080
## 6 5.056433 0.1977679 1.167918 0.9413940 0.6851065 0.8725683
## R_0_S_2 R_0_NGM
## 1 0.003147728 2.138726
## 2 0.019963497 2.108494
## 3 0.004824122 2.350278
## 4 0.023430446 2.270527
## 5 0.002968733 2.492698
## 6 0.008241405 2.507310range(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers$E_0)## [1] 44295.44 71702.50range(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers$z_0)## [1] 9070.932 17695.218range(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers$b_p)## [1] 0.1400547 0.9865475p = ggplot(data = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
aes(x = b_p)) +
geom_histogram() + rahul_theme
p## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
large_b_p_params = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers %>%
filter(b_p >0.90)
p = ggplot(data = large_b_p_params,
aes(x = R_0,
y = R_0_NGM,
color = b_a)) +
geom_point() + scale_color_viridis_c(name = "Asymp \n Trans \n Strength") +
xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number \n (Includes Pre-Symp and Asymp Transmission)") +
rahul_man_figure_theme +
theme_white_background+
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_high_b_p.png")
print(p)
dev.off()## quartz_off_screen
## 2p
small_b_p_params = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers %>%
filter(b_p <0.50)
small_b_p_params## b_a M_0 V_0 K_0 R_0 b_q b_p p_S p_H_cond_S
## 1 0.2413793 5 13 14 8.125304 0.1611435 0.3220870 0.1574198 0.1567379
## 2 0.2758621 5 13 14 7.434055 0.1584788 0.3183492 0.1713446 0.1567255
## 3 0.3448276 5 13 14 6.325641 0.1615416 0.4648878 0.1623848 0.1546026
## 4 0.3448276 5 13 14 6.619633 0.1668452 0.4654342 0.1359111 0.1505722
## 5 0.5517241 5 13 14 5.776447 0.1487695 0.2588775 0.1517047 0.1843260
## 6 0.5862069 5 13 14 6.171486 0.1337064 0.1400547 0.1645830 0.1660020
## 7 0.6896552 5 13 14 4.930404 0.1541091 0.4620370 0.1486280 0.1836424
## phi_E phi_U phi_S h_V gamma N_0 E_0 z_0 C_0
## 1 1.09 1.09 0.2 0.125 79.534133 8e+06 52522.51 12908.19 0
## 2 1.09 1.09 0.2 0.125 3861.955611 8e+06 50636.94 11633.76 0
## 3 1.09 1.09 0.2 0.125 253.068498 8e+06 46679.50 14443.09 0
## 4 1.09 1.09 0.2 0.125 173.699820 8e+06 53744.73 17695.22 0
## 5 1.09 1.09 0.2 0.125 194.903437 8e+06 56052.18 11154.96 0
## 6 1.09 1.09 0.2 0.125 182.122463 8e+06 44295.44 12335.03 0
## 7 1.09 1.09 0.2 0.125 2.663607 8e+06 55541.02 11957.99 0
## social_distancing_start_time quarantine_start_time PCR_sens sigma_M
## 1 17 22 0.9 0.2784369
## 2 17 22 0.9 0.2790559
## 3 17 22 0.9 0.2776869
## 4 17 22 0.9 0.2808862
## 5 17 22 0.9 0.2810877
## 6 17 22 0.9 0.2804963
## 7 17 22 0.9 0.2827579
## beta_w_3 beta_w_2 beta_w_1 beta_w_0 g_0 g_F sigma_epsilon
## 1 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 2 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 3 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 4 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 5 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 6 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 7 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## G_w_y_scaling msg iter_num param_index loglik nfail trace_num
## 1 0.162 mif1 2 121 -628.6123 NA NA
## 2 0.162 mif1 2 138 -628.5431 NA NA
## 3 0.162 mif1 2 168 -629.0010 NA NA
## 4 0.162 mif1 1 180 -629.0174 NA NA
## 5 0.162 mif1 2 271 -628.5576 NA NA
## 6 0.162 mif1 2 288 -629.1750 NA NA
## 7 0.162 mif1 2 343 -629.0228 NA NA
## loglist.se Antibody_Mean_LL Antibody_LL_SE Median_Herd_Immunity
## 1 0.008142806 -24.85219 0.003761241 0.1912882
## 2 0.009305056 -26.00451 0.006860371 0.1764086
## 3 0.010589918 -25.15982 0.003990094 0.1851563
## 4 0.006043226 -24.79345 0.004001890 0.2221081
## 5 0.004688799 -24.69044 0.002171494 0.1986929
## 6 0.007497558 -25.45496 0.004353770 0.1833871
## 7 0.006280933 -24.50344 0.002076568 0.2082881
## combo_num sim_subset_index duration_of_symp_1 duration_of_symp_2
## 1 31 1 5 0.0125732181
## 2 36 1 5 0.0002589362
## 3 43 1 5 0.0039514993
## 4 46 1 5 0.0057570583
## 5 68 1 5 0.0051307458
## 6 70 1 5 0.0054908109
## 7 76 1 5 0.3754308107
## duration_of_symp gamma_total Beta R_0_P R_0_A R_0_S_1
## 1 5.012573 0.1994983 1.6209846 0.4789891 1.648391 1.2758754
## 2 5.000259 0.1999896 1.4867341 0.4342207 1.699297 1.2737192
## 3 5.003951 0.1998421 1.2641291 0.5391543 1.825610 1.0263769
## 4 5.005757 0.1997700 1.3224039 0.5646715 1.970129 0.8986468
## 5 5.005131 0.1997950 1.1541052 0.2741026 2.700750 0.8754160
## 6 5.005491 0.1997806 1.2329432 0.1584215 3.019029 1.0146071
## 7 5.375431 0.1860316 0.9172109 0.3887939 2.692716 0.6816160
## R_0_S_2 R_0_NGM
## 1 0.0027054986 3.405961
## 2 0.0000556244 3.407292
## 3 0.0006857403 3.391827
## 4 0.0008789135 3.434326
## 5 0.0007327260 3.851002
## 6 0.0009292432 4.192987
## 7 0.0417811232 3.804907range(small_b_p_params$b_a)## [1] 0.2413793 0.6896552small_b_p_big_b_a = small_b_p_params %>%
filter(b_a == max(b_a))
small_b_p_small_b_a = small_b_p_params %>%
filter(b_a == min(b_a))
small_b_p_big_b_a$b_p## [1] 0.462037small_b_p_small_b_a$b_p## [1] 0.322087hist(small_b_p_params$b_p)
hist(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers$b_p)
big_b_a_values = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers %>%
filter(b_a >= 0.52)
range(big_b_a_values$R_0_NGM)## [1] 3.153313 4.192987range(big_b_a_values$R_0)## [1] 2.927915 6.171486small_b_a_values =antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers %>%
filter(b_a < 0.52)
range(small_b_a_values$R_0_NGM)## [1] 2.108494 3.434326range(small_b_a_values$R_0)## [1] 3.984487 8.125304range(antibody_top_2_LL_from_b_a_profile_top_2_LL$p_S)## [1] 0.1292286 0.1740948mid_b_p = antibody_top_2_LL_from_b_a_profile_top_2_LL%>%
filter(b_p <0.55) %>%
filter(b_p > 0.45)
mid_b_p_big_b_a = mid_b_p %>%
filter(b_a == max(b_a))
mid_b_p_small_b_a = mid_b_p %>%
filter(b_a < max(b_a)) %>%
filter(b_p == min(b_p))
mid_b_p_big_b_a$R_0_NGM## [1] 3.804907mid_b_p_small_b_a$R_0_NGM## [1] 3.391827mid_b_p_big_b_a$R_0## [1] 4.930404mid_b_p_small_b_a$R_0## [1] 6.325641write.csv(mid_b_p_big_b_a, "../Generated_Data/Profiles/N_12_Model/top_2_LL_data/Man_Table_data/b_a_profile_top_2_LL_via_case_and_antibody_LL_mid_b_p_big_b_a_parameter_combination_for_simulation.csv")
write.csv(mid_b_p_small_b_a, "../Generated_Data/Profiles/N_12_Model/top_2_LL_data/Man_Table_data/b_a_profile_top_2_LL_via_case_and_antibody_LL_mid_b_p_small_b_a_parameter_combination_for_simulation.csv")
low_b_p = antibody_top_2_LL_from_b_a_profile_top_2_LL%>%
filter(b_p < 0.70) %>%
filter(b_p > 0.50)
hist(low_b_p$b_p)
low_b_p$b_a## [1] 0.3448276 0.4827586Pick big b_a parameter combination
big_b_a_MLE = large_b_p_params %>%
filter(b_a == max(b_a)) %>%
filter(b_p == max(b_p))
big_b_a_MLE## b_a M_0 V_0 K_0 R_0 b_q b_p p_S p_H_cond_S
## 1 0.9655172 5 13 14 3.083236 0.163873 0.9865475 0.154388 0.1735144
## phi_E phi_U phi_S h_V gamma N_0 E_0 z_0 C_0
## 1 1.09 1.09 0.2 0.125 11.72791 8e+06 54806.41 11625.01 0
## social_distancing_start_time quarantine_start_time PCR_sens sigma_M
## 1 17 22 0.9 0.27583
## beta_w_3 beta_w_2 beta_w_1 beta_w_0 g_0 g_F sigma_epsilon
## 1 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## G_w_y_scaling msg iter_num param_index loglik nfail trace_num
## 1 0.162 mif1 1 478 -628.7638 NA NA
## loglist.se Antibody_Mean_LL Antibody_LL_SE Median_Herd_Immunity
## 1 0.01206257 -24.50512 0.001865621 0.2021783
## combo_num sim_subset_index duration_of_symp_1 duration_of_symp_2
## 1 94 1 5 0.08526666
## duration_of_symp gamma_total Beta R_0_P R_0_A R_0_S_1
## 1 5.085267 0.1966465 0.6063076 0.5487626 2.475108 0.4680332
## R_0_S_2 R_0_NGM
## 1 0.006596615 3.498501write.csv(big_b_a_MLE, "../Generated_Data/Profiles/N_12_Model/top_2_LL_data/Man_Table_data/b_a_profile_top_2_LL_via_case_and_antibody_LL_big_b_a_parameter_combination_for_simulation.csv")
big_b_a_MLE_comb = big_b_a_MLE %>%
dplyr::select(-msg, -iter_num, -param_index,
-loglik, -nfail, -trace_num, -loglist.se,
-Antibody_Mean_LL, -Antibody_LL_SE,
-Median_Herd_Immunity, -combo_num,
-sim_subset_index, -duration_of_symp_1,
-duration_of_symp_2, -duration_of_symp,
-gamma_total, -Beta,
-R_0_P, -R_0_A, -R_0_S_1, -R_0_S_2,
-R_0_NGM)
mid_b_p_big_b_a_comb = mid_b_p_big_b_a %>%
dplyr::select(-msg, -iter_num, -param_index,
-loglik, -nfail, -trace_num, -loglist.se,
-Antibody_Mean_LL, -Antibody_LL_SE,
-Median_Herd_Immunity, -combo_num,
-sim_subset_index, -duration_of_symp_1,
-duration_of_symp_2, -duration_of_symp,
-gamma_total, -Beta,
-R_0_P, -R_0_A, -R_0_S_1, -R_0_S_2,
-R_0_NGM)Load covaraites and obs data
#Load Observed NYC case data
Observed_data = read.csv(paste0(
"../Generated_Data/observed_data_",
model_name, ".csv"))
head(Observed_data)## Y times obs_prop_positive
## 1 0 1 NA
## 2 0 2 0.00000000
## 3 2 3 0.25000000
## 4 2 4 0.05555556
## 5 7 5 0.15555556
## 6 0 6 0.00000000### Define start date
true_start_date = as.Date("2020-03-01")
t0 = 0
start_of_year = as.Date("2020-01-01")
first_saturday_in_year = as.Date("2020-01-04")
## Compartment/Queue Cohort Numbers
M = 5
V = 13
K = 14
#Declare Csnippets and data
source("Csnippet_nyc_coronavirus_model_N_12.R")
## Load NYC covariate data
covariate_df = read.csv(file =
paste0("../Generated_Data/covariate_data_",
model_name, ".csv"))
### Create covariate table
covar=covariate_table(
time=covariate_df$times,
L_advanced_2_days=covariate_df$L_advanced_2_days,
F_w_y = covariate_df$F_w_y,
L_orig = covariate_df$L_orig,
w = covariate_df$Week,
y = covariate_df$Year,
times="time"
)Simulate from big b_a param combination
##Simulation from ML
sim_data_big_b_a = simulate(nsim = 101,
seed = 12345,
times = Observed_data$times,
t0 = t0,
rprocess = pomp::euler(rproc,delta.t = 1),
params = big_b_a_MLE_comb,
paramnames = paramnames,
statenames = statenames,
obsnames = obsnames,
accumvars = acumvarnames,
rinit = init,
rmeas = rmeas,
partrans = par_trans,
covar = covar,
format = "data.frame")
#head(sim_data)
sim_data_big_a_FOI_data = sim_data_big_b_a %>%
dplyr::select(time, .id, I_P, I_S_1, I_S_2,
I_A, N, beta_t)
sim_data_big_a_FOI_data = sim_data_big_a_FOI_data %>%
mutate(lambda_FOI_S = (beta_t*(I_S_1 + I_S_2))/N,
lambda_FOI_A = (beta_t*big_b_a_MLE$b_a*I_A)/N,
lambda_FOI_P = (beta_t*big_b_a_MLE$b_p*I_P)/N) %>%
mutate(lambda_FOI_total =
lambda_FOI_S + lambda_FOI_A + lambda_FOI_P)
median_FOI_data_big_b_a = sim_data_big_a_FOI_data %>%
group_by(time) %>%
summarize(Symptomatic = median(lambda_FOI_S),
Asymptomatic = median(lambda_FOI_A),
Presymptomatic = median(lambda_FOI_P)) %>%
as.data.frame()
median_FOI_data_big_b_a_melt = median_FOI_data_big_b_a %>%
melt(id.vars= c("time")) %>%
dplyr::select(time, FOI = value,
Infection_Status = variable)
#head(median_FOI_data)
p = ggplot(data = median_FOI_data_big_b_a_melt,
aes(x = time,
y = FOI,
fill = Infection_Status)) +
geom_area() + rahul_man_figure_theme +
theme_white_background +
xlab("Days since March 1, 2020") +
ylab("Force of Infection") +
theme(legend.position = c(0.70,0.75)) +
scale_fill_discrete(name = "Infection Status")+
theme(axis.title.x = element_text(face = "plain", size = 24),
axis.title.y = element_text(face = "plain", size = 24))+
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_large_b_a_parameter_combination.png")
print(p)
dev.off()## quartz_off_screen
## 2jpeg("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_large_b_a_parameter_combination.jpeg",
quality = 100)
print(p)
dev.off()## quartz_off_screen
## 2png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_large_b_a_parameter_combination.png")
print(p)
dev.off()## quartz_off_screen
## 2jpeg("../Figures/Profiles/N_12_Model/Man_Figs/contribution_to_FOI_for_large_b_a_parameter_combination.jpeg",
quality = 100)
print(p)
dev.off()## quartz_off_screen
## 2png("../Figures/Profiles/N_12_Model/Man_Figs/contribution_to_FOI_for_large_b_a_parameter_combination.png")
print(p)
dev.off()## quartz_off_screen
## 2Simulate from mid b_p big b_a param combination
##Simulation from ML
sim_data_mid_b_p_big_b_a = simulate(nsim = 101,
seed = 12345,
times = Observed_data$times,
t0 = t0,
rprocess = pomp::euler(rproc,delta.t = 1),
params = mid_b_p_big_b_a_comb,
paramnames = paramnames,
statenames = statenames,
obsnames = obsnames,
accumvars = acumvarnames,
rinit = init,
rmeas = rmeas,
partrans = par_trans,
covar = covar,
format = "data.frame")
#head(sim_data)
sim_data_mid_b_p_big_a_FOI_data = sim_data_mid_b_p_big_b_a %>%
dplyr::select(time, .id, I_P, I_S_1, I_S_2,
I_A, N, beta_t)
sim_data_mid_b_p_big_a_FOI_data = sim_data_mid_b_p_big_a_FOI_data %>%
mutate(lambda_FOI_S = (beta_t*(I_S_1 + I_S_2))/N,
lambda_FOI_A = (beta_t*mid_b_p_big_b_a_comb$b_a*I_A)/N,
lambda_FOI_P = (beta_t*mid_b_p_big_b_a_comb$b_p*I_P)/N) %>%
mutate(lambda_FOI_total =
lambda_FOI_S + lambda_FOI_A + lambda_FOI_P)
median_FOI_data_mid_b_p_big_b_a = sim_data_mid_b_p_big_a_FOI_data %>%
group_by(time) %>%
summarize(Symptomatic = median(lambda_FOI_S),
Asymptomatic = median(lambda_FOI_A),
Presymptomatic = median(lambda_FOI_P)) %>%
as.data.frame()
median_FOI_data_mid_b_p_big_b_a_melt = median_FOI_data_mid_b_p_big_b_a %>%
melt(id.vars= c("time")) %>%
dplyr::select(time, FOI = value,
Infection_Status = variable)
#head(median_FOI_data)
p = ggplot(data = median_FOI_data_mid_b_p_big_b_a_melt,
aes(x = time,
y = FOI,
fill = Infection_Status)) +
geom_area() + rahul_man_figure_theme +
theme_white_background +
xlab("Days since March 1, 2020") +
ylab("Force of Infection") +
theme(legend.position = c(0.70,0.75)) +
scale_fill_discrete(name = "Infection Status")+
theme(axis.title.x = element_text(face = "plain", size = 24),
axis.title.y = element_text(face = "plain", size = 24))+
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_mid_b_p_large_b_a_parameter_combination.png")
print(p)
dev.off()## quartz_off_screen
## 2jpeg("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_mid_b_p_large_b_a_parameter_combination.jpeg",
quality = 100)
print(p)
dev.off()## quartz_off_screen
## 2png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_mid_b_p_large_b_a_parameter_combination.png")
print(p)
dev.off()## quartz_off_screen
## 2jpeg("../Figures/Profiles/N_12_Model/Man_Figs/contribution_to_FOI_for_mid_b_p_large_b_a_parameter_combination.jpeg",
quality = 100)
print(p)
dev.off()## quartz_off_screen
## 2png("../Figures/Profiles/N_12_Model/Man_Figs/contribution_to_FOI_for_mid_b_p_large_b_a_parameter_combination.png")
print(p)
dev.off()## quartz_off_screen
## 2Identify small b_a param
small_b_a_MLE = large_b_p_params %>%
filter(b_a > 0) %>%
filter(b_a == min(b_a)) %>%
filter(b_p == max(b_p))
small_b_a_MLE## b_a M_0 V_0 K_0 R_0 b_q b_p p_S p_H_cond_S
## 1 0.06896552 5 13 14 6.098822 0.2343746 0.940249 0.148861 0.1569547
## phi_E phi_U phi_S h_V gamma N_0 E_0 z_0 C_0
## 1 1.09 1.09 0.2 0.125 6.333218 8e+06 63566.34 13443.03 0
## social_distancing_start_time quarantine_start_time PCR_sens sigma_M
## 1 17 22 0.9 0.2709327
## beta_w_3 beta_w_2 beta_w_1 beta_w_0 g_0 g_F sigma_epsilon
## 1 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## G_w_y_scaling msg iter_num param_index loglik nfail trace_num
## 1 0.162 mif1 2 44 -627.43 NA NA
## loglist.se Antibody_Mean_LL Antibody_LL_SE Median_Herd_Immunity
## 1 0.008839892 -24.41689 0.002390176 0.1998329
## combo_num sim_subset_index duration_of_symp_1 duration_of_symp_2
## 1 13 1 5 0.1578976
## duration_of_symp gamma_total Beta R_0_P R_0_A R_0_S_1
## 1 5.157898 0.1938774 1.182424 1.019975 0.347037 0.880084
## R_0_S_2 R_0_NGM
## 1 0.02343045 2.270527write.csv(small_b_a_MLE, "../Generated_Data/Profiles/N_12_Model/top_2_LL_data/Man_Table_data/b_a_profile_top_2_LL_via_case_and_antibody_LL_small_b_a_parameter_combination_for_simulation.csv")
small_b_a_MLE_comb = small_b_a_MLE %>%
dplyr::select(-msg, -iter_num, -param_index,
-loglik, -nfail, -trace_num, -loglist.se,
-Antibody_Mean_LL, -Antibody_LL_SE,
-Median_Herd_Immunity, -combo_num,
-sim_subset_index, -duration_of_symp_1,
-duration_of_symp_2, -duration_of_symp,
-gamma_total, -Beta,
-R_0_P, -R_0_A, -R_0_S_1, -R_0_S_2,
-R_0_NGM)Identify mid b_p small b_a param
mid_b_p_small_b_a_comb = mid_b_p_small_b_a %>%
dplyr::select(-msg, -iter_num, -param_index,
-loglik, -nfail, -trace_num, -loglist.se,
-Antibody_Mean_LL, -Antibody_LL_SE,
-Median_Herd_Immunity, -combo_num,
-sim_subset_index, -duration_of_symp_1,
-duration_of_symp_2, -duration_of_symp,
-gamma_total, -Beta,
-R_0_P, -R_0_A, -R_0_S_1, -R_0_S_2,
-R_0_NGM)Simulate from small b_a param combination
##Simulation from ML
sim_data_small_b_a = simulate(nsim = 101,
seed = 12345,
times = Observed_data$times,
t0 = t0,
rprocess = pomp::euler(rproc,delta.t = 1),
params = small_b_a_MLE_comb,
paramnames = paramnames,
statenames = statenames,
obsnames = obsnames,
accumvars = acumvarnames,
rinit = init,
rmeas = rmeas,
partrans = par_trans,
covar = covar,
format = "data.frame")
#head(sim_data)
sim_data_small_a_FOI_data = sim_data_small_b_a %>%
dplyr::select(time, .id, I_P, I_S_1, I_S_2,
I_A, N, beta_t)
sim_data_small_a_FOI_data = sim_data_small_a_FOI_data %>%
mutate(lambda_FOI_S = (beta_t*(I_S_1 + I_S_2))/N,
lambda_FOI_A = (beta_t*small_b_a_MLE$b_a*I_A)/N,
lambda_FOI_P = (beta_t*small_b_a_MLE$b_p*I_P)/N) %>%
mutate(lambda_FOI_total =
lambda_FOI_S + lambda_FOI_A + lambda_FOI_P)
median_FOI_data_small_b_a = sim_data_small_a_FOI_data %>%
group_by(time) %>%
summarize(Symptomatic = median(lambda_FOI_S),
Asymptomatic = median(lambda_FOI_A),
Presymptomatic = median(lambda_FOI_P)) %>%
as.data.frame()
median_FOI_data_small_b_a_melt = median_FOI_data_small_b_a %>%
melt(id.vars= c("time")) %>%
dplyr::select(time, FOI = value,
Infection_Status = variable)
#head(median_FOI_data)
p = ggplot(data = median_FOI_data_small_b_a_melt,
aes(x = time,
y = FOI,
fill = Infection_Status)) +
geom_area() + rahul_man_figure_theme +
theme_white_background +
xlab("Days since March 1, 2020") +
ylab("Force of Infection") +
theme(legend.position = c(0.70,0.75)) +
scale_fill_discrete(name = "Infection Status")+
theme(axis.title.x = element_text(face = "plain", size = 24),
axis.title.y = element_text(face = "plain", size = 24))+
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_small_b_a_parameter_combination.png")
print(p)
dev.off()## quartz_off_screen
## 2jpeg("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_small_b_a_parameter_combination.jpeg",
quality = 100)
print(p)
dev.off()## quartz_off_screen
## 2jpeg("../Figures/Profiles/N_12_Model/Man_Figs/contribution_to_FOI_for_small_b_a_parameter_combination.jpeg",
quality = 100)
print(p)
dev.off()## quartz_off_screen
## 2png("../Figures/Profiles/N_12_Model/Man_Figs/contribution_to_FOI_for_small_b_a_parameter_combination.png")
print(p)
dev.off()## quartz_off_screen
## 2Simulate from small b_a param combination
##Simulation from ML
sim_data_mid_b_p_small_b_a = simulate(nsim = 101,
seed = 12345,
times = Observed_data$times,
t0 = t0,
rprocess = pomp::euler(rproc,delta.t = 1),
params = mid_b_p_small_b_a_comb,
paramnames = paramnames,
statenames = statenames,
obsnames = obsnames,
accumvars = acumvarnames,
rinit = init,
rmeas = rmeas,
partrans = par_trans,
covar = covar,
format = "data.frame")
#head(sim_data)
sim_data_mid_b_p_small_a_FOI_data = sim_data_mid_b_p_small_b_a %>%
dplyr::select(time, .id, I_P, I_S_1, I_S_2,
I_A, N, beta_t)
sim_data_mid_b_p_small_a_FOI_data = sim_data_mid_b_p_small_a_FOI_data %>%
mutate(lambda_FOI_S = (beta_t*(I_S_1 + I_S_2))/N,
lambda_FOI_A = (beta_t*mid_b_p_small_b_a_comb$b_a*I_A)/N,
lambda_FOI_P = (beta_t*mid_b_p_small_b_a_comb$b_p*I_P)/N) %>%
mutate(lambda_FOI_total =
lambda_FOI_S + lambda_FOI_A + lambda_FOI_P)
median_FOI_data_mid_b_p_small_b_a = sim_data_mid_b_p_small_a_FOI_data %>%
group_by(time) %>%
summarize(Symptomatic = median(lambda_FOI_S),
Asymptomatic = median(lambda_FOI_A),
Presymptomatic = median(lambda_FOI_P)) %>%
as.data.frame()
median_FOI_data_mid_b_p_small_b_a_melt = median_FOI_data_mid_b_p_small_b_a %>%
melt(id.vars= c("time")) %>%
dplyr::select(time, FOI = value,
Infection_Status = variable)
#head(median_FOI_data)
p = ggplot(data = median_FOI_data_mid_b_p_small_b_a_melt,
aes(x = time,
y = FOI,
fill = Infection_Status)) +
geom_area() + rahul_man_figure_theme +
theme_white_background +
xlab("Days since March 1, 2020") +
ylab("Force of Infection") +
theme(legend.position = c(0.70,0.75)) +
scale_fill_discrete(name = "Infection Status")+
theme(axis.title.x = element_text(face = "plain", size = 24),
axis.title.y = element_text(face = "plain", size = 24))+
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_mid_b_p_small_b_a_parameter_combination.png")
print(p)
dev.off()## quartz_off_screen
## 2jpeg("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/contribution_to_FOI_for_mid_b_p_small_b_a_parameter_combination.jpeg",
quality = 100)
print(p)
dev.off()## quartz_off_screen
## 2jpeg("../Figures/Profiles/N_12_Model/Man_Figs/contribution_to_FOI_for_mid_b_p_small_b_a_parameter_combination.jpeg",
quality = 100)
print(p)
dev.off()## quartz_off_screen
## 2png("../Figures/Profiles/N_12_Model/Man_Figs/contribution_to_FOI_for_mid_b_p_small_b_a_parameter_combination.png")
print(p)
dev.off()## quartz_off_screen
## 2outlier = antibody_top_2_LL_from_b_a_profile_top_2_LL %>%
filter(R_0 > 15)
outlier## b_a M_0 V_0 K_0 R_0 b_q b_p p_S
## 1 0.06896552 5 13 14 17.77109 0.133977 0.008385475 0.1630488
## p_H_cond_S phi_E phi_U phi_S h_V gamma N_0 E_0 z_0 C_0
## 1 0.1729462 1.09 1.09 0.2 0.125 960.6433 8e+06 54219.45 11091.14 0
## social_distancing_start_time quarantine_start_time PCR_sens sigma_M
## 1 17 22 0.9 0.2807117
## beta_w_3 beta_w_2 beta_w_1 beta_w_0 g_0 g_F sigma_epsilon
## 1 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## G_w_y_scaling msg iter_num param_index loglik nfail trace_num
## 1 0.162 mif1 1 37 -629.0826 NA NA
## loglist.se Antibody_Mean_LL Antibody_LL_SE Median_Herd_Immunity
## 1 0.004121517 -25.32007 0.007746291 0.1847238
## combo_num sim_subset_index duration_of_symp_1 duration_of_symp_2
## 1 10 1 5 0.001040969
## duration_of_symp gamma_total Beta R_0_P R_0_A R_0_S_1
## 1 5.001041 0.1999584 3.553478 0.02733725 1.025547 2.896952
## R_0_S_2 R_0_NGM
## 1 0.0004988189 3.950335Plot b_a * b_p
product_data =
antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers %>%
mutate(b_a_times_b_p = b_a*b_p,
b_a_plus_b_p = b_a + b_p)
p = ggplot(data = product_data,
aes(x = R_0,
y = R_0_NGM,
color = b_a_times_b_p)) + rahul_man_figure_theme +
scale_color_viridis_c() +
geom_point(size = 5) +
scale_x_continuous(breaks=c(seq(2,9))) +
scale_y_continuous(breaks=seq(2,5,1)) +
coord_cartesian(expand = FALSE,
xlim = c(1.75, 9), ylim = c(1.75, 5.25)) +
xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21)) +
theme(axis.title.x = element_text(face = "plain",
size = 24),
axis.title.y = element_text(face = "plain",
size = 24))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_times_b_p.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = product_data,
aes(x = R_0,
y = R_0_NGM,
color = b_a_times_b_p)) + rahul_man_figure_theme +
scale_color_viridis_c() +
geom_point(size = 5) +
scale_x_continuous(breaks=c(seq(2,9))) +
scale_y_continuous(breaks=seq(2,5,1)) +
coord_cartesian(expand = FALSE,
xlim = c(1.75, 9), ylim = c(1.75, 5.25)) +
xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21)) +
theme(axis.title.x = element_text(face = "plain",
size = 24),
axis.title.y = element_text(face = "plain",
size = 24))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_times_b_p_circle_mid_b_p.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = product_data,
aes(x = R_0,
y = R_0_NGM,
color = b_a_plus_b_p)) + rahul_man_figure_theme +
scale_color_viridis_c() +
geom_point() +
xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21)) +
theme(axis.title.x = element_text(face = "plain",
size = 24),
axis.title.y = element_text(face = "plain",
size = 24))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_plus_b_p.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = product_data,
aes(x = R_0,
y = R_0_NGM,
color = b_a_plus_b_p)) + rahul_man_figure_theme +
scale_color_viridis_c() +
geom_point() +
xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_plus_b_p_circle_mid_b_p.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = product_data,
aes(x = R_0,
y = R_0_NGM,
color = b_a)) +
rahul_man_figure_theme +
scale_color_viridis_c(name = expression(b[a])) +
geom_point(size = 5) +
scale_x_continuous(breaks=c(seq(2,9))) +
scale_y_continuous(breaks=seq(2,5,1)) +
coord_cartesian(expand = FALSE,
xlim = c(1.75, 9), ylim = c(1.75, 5.25)) +
theme_white_background +
xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
theme(axis.title.x = element_text(face = "plain",
size = 24),
axis.title.y = element_text(face = "plain",
size = 24))+
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21))#+ #+
# theme(legend.position = c(0.75,0.75))
# labs(y = expression(R["0 NGM"]),
# x = expression(R[0]))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = product_data,
aes(x = R_0,
y = R_0_NGM,
color = b_a)) +
rahul_man_figure_theme +
scale_color_viridis_c(name = expression(b[a])) +
geom_point(size = 5) +
scale_x_continuous(breaks=c(seq(2,9))) +
scale_y_continuous(breaks=seq(2,5,1)) +
coord_cartesian(expand = FALSE,
xlim = c(1.75, 9), ylim = c(1.75, 5.25)) +
theme_white_background +
xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21))+ #+
# theme(legend.position = c(0.75,0.75))
labs(y = expression(R[0["NGM"]]),
x = expression(R[0]))+
theme(axis.title.x = element_text(face = "plain",
size = 32),
axis.title.y = element_text(face = "plain",
size = 32))
p
#expression(R[0["NGM"]])
png("../Figures/Profiles/N_12_Model/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_symb_plot.png")
print(p)
dev.off()## quartz_off_screen
## 2write.csv(product_data,
"../Generated_Data/Profiles/N_12_Model/top_2_LL_data/product_data_for_Fig_3C.csv", row.names = FALSE)
p = ggplot(data = product_data,
aes(x = R_0,
y = R_0_NGM,
color = b_a)) +
rahul_man_figure_theme +
scale_color_viridis_c(name = expression(b[a])) +
geom_point(size = 5) +
scale_x_continuous(breaks=c(seq(2,9))) +
scale_y_continuous(breaks=seq(2,5,1)) +
coord_cartesian(expand = FALSE,
xlim = c(1.75, 9), ylim = c(1.75, 5.25)) +
theme_white_background +
xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
theme(axis.title.x = element_text(face = "plain", size = 24),
axis.title.y = element_text(face = "plain", size = 24))+
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21)) #+
# theme(legend.position = c(0.75,0.75))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_circle_mid_b_p.png")
print(p)
dev.off()## quartz_off_screen
## 2png("../Figures/Profiles/N_12_Model/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_circle_mid_b_p.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = product_data,
aes(x = R_0,
y = R_0_NGM,
color = b_p)) +
rahul_man_figure_theme +
scale_color_viridis_c() +
geom_point() +
xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_p.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = product_data,
aes(x = R_0,
y = R_0_NGM,
color = b_p)) +
rahul_man_figure_theme +
scale_color_viridis_c() +
geom_point() +
xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_p_circle_mid_b_p.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = product_data,
aes(x = R_0,
y = R_0_NGM,
color = b_a*b_p)) +
rahul_man_figure_theme +
scale_color_viridis_c(name = expression(paste(b[a], "*", b[p]))) +
geom_point(size = 5) +
scale_x_continuous(breaks=c(seq(2,9))) +
scale_y_continuous(breaks=seq(2,5,1)) +
coord_cartesian(expand = FALSE,
xlim = c(1.75, 9), ylim = c(1.75, 5.25)) +
xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
theme_white_background +
theme(axis.title.x = element_text(face = "plain", size = 24),
axis.title.y = element_text(face = "plain", size = 24)) +
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21))
p
png("../Figures/Profiles/N_12_Model/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_times_b_p.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = product_data,
aes(x = R_0,
y = R_0_NGM,
color = b_a*b_p)) +
rahul_man_figure_theme +
scale_color_viridis_c(name = expression(paste(b[a], "*", b[p]))) +
geom_point(size = 5) +
scale_x_continuous(breaks=c(seq(2,9))) +
scale_y_continuous(breaks=seq(2,5,1)) +
coord_cartesian(expand = FALSE,
xlim = c(1.75, 9), ylim = c(1.75, 5.25)) +
xlab("Secondary cases from each \n primary symptomatic case") + ylab("Overall Reproductive Number") +
theme_white_background +
theme(axis.title.x = element_text(face = "plain", size = 24),
axis.title.y = element_text(face = "plain", size = 24))+
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21)) #+
# theme(legend.position = c(0.75, 0.75))
p
png("../Figures/Profiles/N_12_Model/Man_Figs/b_a_profile_R_0_vs_R_0_NGM_antibody_LL_top_no_outliers_color_b_a_times_b_p_circle_mid_b_p.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = product_data,
aes(x = b_a*b_p,
y = R_0)) +
rahul_man_figure_theme +
geom_point(size= 5) +
ylab("Secondary cases from each \n primary symptomatic case") +
xlab(expression(paste(b[a],'*',b[p]))) +
theme_white_background+
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_b_a_times_b_p_vs_R_0_antibody_LL_top_no_outliers.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = product_data,
aes(x = b_a*b_p,
y = R_0)) +
rahul_man_figure_theme +
geom_point(size= 5) +
ylab("Secondary cases from each \n primary symptomatic case") +
xlab(expression(paste(b[a],'*',b[p]))) +
theme_white_background+
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_b_a_times_b_p_vs_R_0_antibody_LL_top_no_outliers_circle_mid_b_p.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = product_data,
aes(x = b_a+b_p,
y = R_0)) +
rahul_man_figure_theme +
geom_point() +
ylab("Secondary cases from each \n primary symptomatic case") +
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_b_a_plus_b_p_vs_R_0_antibody_LL_top_no_outliers.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = product_data,
aes(x = b_a+b_p,
y = R_0)) +
rahul_man_figure_theme +
geom_point() +
ylab("Secondary cases from each \n primary symptomatic case")+
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) +
theme(axis.text.x = element_text(size=21)) +
theme(axis.text.y = element_text(size=21))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/top_2_LL_via_antibody_comp_plots/Man_Figs/b_a_profile_b_a_plus_b_p_vs_R_0_antibody_LL_top_no_outliers_cirlce_mid_b_p.png")
print(p)
dev.off()## quartz_off_screen
## 2head(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers)## b_a M_0 V_0 K_0 R_0 b_q b_p p_S p_H_cond_S
## 1 0.00000000 5 13 14 7.061999 0.2401933 0.7957984 0.1574275 0.1749811
## 2 0.00000000 5 13 14 6.215237 0.2395227 0.9233632 0.1740948 0.1532892
## 3 0.06896552 5 13 14 5.984521 0.2204782 0.8913949 0.1735879 0.1757053
## 4 0.06896552 5 13 14 6.098822 0.2343746 0.9402490 0.1488610 0.1569547
## 5 0.10344828 5 13 14 6.845525 0.2237219 0.7932278 0.1296082 0.1818042
## 6 0.13793103 5 13 14 5.905499 0.2168981 0.8785885 0.1494229 0.1631618
## phi_E phi_U phi_S h_V gamma N_0 E_0 z_0 C_0
## 1 1.09 1.09 0.2 0.125 58.078038 8e+06 65553.63 10267.580 0
## 2 1.09 1.09 0.2 0.125 8.978508 8e+06 59637.11 11293.462 0
## 3 1.09 1.09 0.2 0.125 35.301204 8e+06 56854.62 9411.456 0
## 4 1.09 1.09 0.2 0.125 6.333218 8e+06 63566.34 13443.034 0
## 5 1.09 1.09 0.2 0.125 48.705248 8e+06 71702.50 12381.941 0
## 6 1.09 1.09 0.2 0.125 17.720241 8e+06 63878.93 12506.505 0
## social_distancing_start_time quarantine_start_time PCR_sens sigma_M
## 1 17 22 0.9 0.2764235
## 2 17 22 0.9 0.2750787
## 3 17 22 0.9 0.2768054
## 4 17 22 0.9 0.2709327
## 5 17 22 0.9 0.2752184
## 6 17 22 0.9 0.2751406
## beta_w_3 beta_w_2 beta_w_1 beta_w_0 g_0 g_F sigma_epsilon
## 1 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 2 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 3 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 4 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 5 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 6 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## G_w_y_scaling msg iter_num param_index loglik nfail trace_num
## 1 0.162 mif1 1 10 -627.3202 NA NA
## 2 0.162 mif1 1 11 -627.2093 NA NA
## 3 0.162 mif1 2 43 -628.1031 NA NA
## 4 0.162 mif1 2 44 -627.4300 NA NA
## 5 0.162 mif1 1 60 -628.1314 NA NA
## 6 0.162 mif1 2 76 -628.5305 NA NA
## loglist.se Antibody_Mean_LL Antibody_LL_SE Median_Herd_Immunity
## 1 0.008206760 -24.75616 0.004776843 0.1933564
## 2 0.010463251 -26.12591 0.007764530 0.1749998
## 3 0.007808338 -26.08064 0.008118650 0.1775874
## 4 0.008839892 -24.41689 0.002390176 0.1998329
## 5 0.010838884 -25.56156 0.006393602 0.2380878
## 6 0.008314098 -24.34239 0.002107043 0.2090317
## combo_num sim_subset_index duration_of_symp_1 duration_of_symp_2
## 1 3 1 5 0.01721821
## 2 4 1 5 0.11137708
## 3 11 1 5 0.02832765
## 4 13 1 5 0.15789761
## 5 15 1 5 0.02053167
## 6 19 1 5 0.05643264
## duration_of_symp gamma_total Beta R_0_P R_0_A R_0_S_1
## 1 5.017218 0.1993136 1.407553 1.0276405 0.0000000 1.1079374
## 2 5.111377 0.1956420 1.215961 1.0300678 0.0000000 1.0584630
## 3 5.028328 0.1988733 1.190161 0.9733061 0.3391599 1.0329881
## 4 5.157898 0.1938774 1.182424 1.0199752 0.3470370 0.8800840
## 5 5.020532 0.1991821 1.363506 0.9922669 0.6138540 0.8836080
## 6 5.056433 0.1977679 1.167918 0.9413940 0.6851065 0.8725683
## R_0_S_2 R_0_NGM
## 1 0.003147728 2.138726
## 2 0.019963497 2.108494
## 3 0.004824122 2.350278
## 4 0.023430446 2.270527
## 5 0.002968733 2.492698
## 6 0.008241405 2.507310p = ggplot(data = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
aes(x = p_H_cond_S, y = b_a)) +
geom_point() + rahul_man_figure_theme
p
p = ggplot(data = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
aes(x = p_S, y = b_a)) +
geom_point() + rahul_man_figure_theme
p
head(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers)## b_a M_0 V_0 K_0 R_0 b_q b_p p_S p_H_cond_S
## 1 0.00000000 5 13 14 7.061999 0.2401933 0.7957984 0.1574275 0.1749811
## 2 0.00000000 5 13 14 6.215237 0.2395227 0.9233632 0.1740948 0.1532892
## 3 0.06896552 5 13 14 5.984521 0.2204782 0.8913949 0.1735879 0.1757053
## 4 0.06896552 5 13 14 6.098822 0.2343746 0.9402490 0.1488610 0.1569547
## 5 0.10344828 5 13 14 6.845525 0.2237219 0.7932278 0.1296082 0.1818042
## 6 0.13793103 5 13 14 5.905499 0.2168981 0.8785885 0.1494229 0.1631618
## phi_E phi_U phi_S h_V gamma N_0 E_0 z_0 C_0
## 1 1.09 1.09 0.2 0.125 58.078038 8e+06 65553.63 10267.580 0
## 2 1.09 1.09 0.2 0.125 8.978508 8e+06 59637.11 11293.462 0
## 3 1.09 1.09 0.2 0.125 35.301204 8e+06 56854.62 9411.456 0
## 4 1.09 1.09 0.2 0.125 6.333218 8e+06 63566.34 13443.034 0
## 5 1.09 1.09 0.2 0.125 48.705248 8e+06 71702.50 12381.941 0
## 6 1.09 1.09 0.2 0.125 17.720241 8e+06 63878.93 12506.505 0
## social_distancing_start_time quarantine_start_time PCR_sens sigma_M
## 1 17 22 0.9 0.2764235
## 2 17 22 0.9 0.2750787
## 3 17 22 0.9 0.2768054
## 4 17 22 0.9 0.2709327
## 5 17 22 0.9 0.2752184
## 6 17 22 0.9 0.2751406
## beta_w_3 beta_w_2 beta_w_1 beta_w_0 g_0 g_F sigma_epsilon
## 1 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 2 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 3 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 4 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 5 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## 6 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005 109.1121
## G_w_y_scaling msg iter_num param_index loglik nfail trace_num
## 1 0.162 mif1 1 10 -627.3202 NA NA
## 2 0.162 mif1 1 11 -627.2093 NA NA
## 3 0.162 mif1 2 43 -628.1031 NA NA
## 4 0.162 mif1 2 44 -627.4300 NA NA
## 5 0.162 mif1 1 60 -628.1314 NA NA
## 6 0.162 mif1 2 76 -628.5305 NA NA
## loglist.se Antibody_Mean_LL Antibody_LL_SE Median_Herd_Immunity
## 1 0.008206760 -24.75616 0.004776843 0.1933564
## 2 0.010463251 -26.12591 0.007764530 0.1749998
## 3 0.007808338 -26.08064 0.008118650 0.1775874
## 4 0.008839892 -24.41689 0.002390176 0.1998329
## 5 0.010838884 -25.56156 0.006393602 0.2380878
## 6 0.008314098 -24.34239 0.002107043 0.2090317
## combo_num sim_subset_index duration_of_symp_1 duration_of_symp_2
## 1 3 1 5 0.01721821
## 2 4 1 5 0.11137708
## 3 11 1 5 0.02832765
## 4 13 1 5 0.15789761
## 5 15 1 5 0.02053167
## 6 19 1 5 0.05643264
## duration_of_symp gamma_total Beta R_0_P R_0_A R_0_S_1
## 1 5.017218 0.1993136 1.407553 1.0276405 0.0000000 1.1079374
## 2 5.111377 0.1956420 1.215961 1.0300678 0.0000000 1.0584630
## 3 5.028328 0.1988733 1.190161 0.9733061 0.3391599 1.0329881
## 4 5.157898 0.1938774 1.182424 1.0199752 0.3470370 0.8800840
## 5 5.020532 0.1991821 1.363506 0.9922669 0.6138540 0.8836080
## 6 5.056433 0.1977679 1.167918 0.9413940 0.6851065 0.8725683
## R_0_S_2 R_0_NGM
## 1 0.003147728 2.138726
## 2 0.019963497 2.108494
## 3 0.004824122 2.350278
## 4 0.023430446 2.270527
## 5 0.002968733 2.492698
## 6 0.008241405 2.507310p = ggplot(data = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
aes(x = p_S)) +
geom_histogram() +
rahul_man_figure_theme
p## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
range(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers$p_S)## [1] 0.1292286 0.1740948Averaging over time
FOI_df = data.frame(matrix(nrow = 0, ncol =
ncol( antibody_top_2_LL_from_b_a_profile_top_2_LL) + 7))
colnames(FOI_df) = c(colnames(
antibody_top_2_LL_from_b_a_profile_top_2_LL
), "Mean_Symp_prop", "Mean_Asymp_prop", "Mean_Presymp_Prop", "Mean_Symp", "Mean_Asymp", "Mean_Presymp", "Average_Total_FOI")
for(param_index in seq(1:nrow(antibody_top_2_LL_from_b_a_profile_top_2_LL))){
print(param_index)
param_combo =
antibody_top_2_LL_from_b_a_profile_top_2_LL[param_index,]
param_combo_for_sim = param_combo %>%
dplyr::select(-msg, -iter_num, -param_index,
-loglik, -nfail, -trace_num, -loglist.se,
-Antibody_Mean_LL, -Antibody_LL_SE,
-Median_Herd_Immunity, -combo_num,
-sim_subset_index, -duration_of_symp_1,
-duration_of_symp_2, -duration_of_symp,
-gamma_total, -Beta,
-R_0_P, -R_0_A, -R_0_S_1, -R_0_S_2,
-R_0_NGM)
sim_data_big_single_param = simulate(nsim = 101,
seed = 12345,
times = Observed_data$times,
t0 = t0,
rprocess = pomp::euler(rproc,delta.t = 1),
params = param_combo_for_sim,
paramnames = paramnames,
statenames = statenames,
obsnames = obsnames,
accumvars = acumvarnames,
rinit = init,
rmeas = rmeas,
partrans = par_trans,
covar = covar,
format = "data.frame")
#head(sim_data)
sim_data_single_param_FOI_data = sim_data_big_single_param %>%
dplyr::select(time, .id, I_P, I_S_1, I_S_2,
I_A, N, beta_t)
sim_data_single_param_FOI_data = sim_data_single_param_FOI_data %>%
mutate(lambda_FOI_S = (beta_t*(I_S_1 + I_S_2))/N,
lambda_FOI_A = (beta_t*param_combo_for_sim$b_a*I_A)/N,
lambda_FOI_P = (beta_t*param_combo_for_sim$b_p*I_P)/N) %>%
mutate(lambda_FOI_total =
lambda_FOI_S + lambda_FOI_A + lambda_FOI_P)
median_FOI_data_single_param = sim_data_single_param_FOI_data %>%
group_by(time) %>%
summarize(Symptomatic = median(lambda_FOI_S),
Asymptomatic = median(lambda_FOI_A),
Presymptomatic = median(lambda_FOI_P)) %>%
as.data.frame() %>%
mutate(Total = Symptomatic + Asymptomatic + Presymptomatic) %>%
mutate(Symp_Prop = Symptomatic/Total,
Asymp_Prop = Asymptomatic/Total,
Presymp_Prop = Presymptomatic/Total)
single_param_avg_FOI = median_FOI_data_single_param %>%
summarize(Mean_Symp_prop = mean(Symp_Prop),
Mean_Asymp_Prop = mean(Asymp_Prop),
Mean_Presymp_Prop = mean(Presymp_Prop),
Mean_Symp = mean(Symptomatic),
Mean_Asymp = mean(Asymptomatic),
Mean_Presymp = mean(Presymptomatic),
Average_Total_FOI = mean(Total))
single_param_FOI_df = cbind(param_combo, single_param_avg_FOI)
FOI_df = rbind(FOI_df,single_param_FOI_df)
}## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21
## [1] 22
## [1] 23
## [1] 24
## [1] 25
## [1] 26
## [1] 27
## [1] 28
## [1] 29
## [1] 30
## [1] 31
## [1] 32
## [1] 33
## [1] 34
## [1] 35
## [1] 36p = ggplot(data = FOI_df,
aes(x = b_a, y = b_p,
color = Mean_Asymp_Prop)) +
geom_point() + scale_color_viridis_c() +
rahul_man_figure_theme + theme_white_background
p
p = ggplot(data = FOI_df,
aes(x = b_a, y = Mean_Symp)) +
geom_point() +
rahul_man_figure_theme + theme_white_background
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Mean_Symp_FOI.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = FOI_df,
aes(x = b_a, y = Mean_Presymp)) +
geom_point() +
rahul_man_figure_theme + theme_white_background
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Mean_Presymp_FOI.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = FOI_df,
aes(x = b_a, y = Mean_Asymp)) +
geom_point() +
rahul_man_figure_theme + theme_white_background
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Mean_Asymp_FOI.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = FOI_df,
aes(x = b_a, y = Average_Total_FOI)) +
geom_point() +
rahul_man_figure_theme + theme_white_background
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Average_Total_FOI.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = FOI_df,
aes(x = b_a, y = b_p,
color = Mean_Symp_prop)) +
geom_point() + scale_color_viridis_c() +
rahul_man_figure_theme + theme_white_background
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_b_p_color_Mean_symp_prop_of_FOI.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = FOI_df,
aes(x = b_a, y = b_p,
color = (R_0_S_1 + R_0_S_2)/R_0_NGM)) +
geom_point() + scale_color_viridis_c(name = "Prop Symp R_0") +
rahul_man_figure_theme + theme_white_background
p
p = ggplot(data = FOI_df,
aes(x = b_a, y = Mean_Symp_prop)) +
geom_point() +
rahul_man_figure_theme + theme_white_background
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Mean_symp_prop_of_FOI.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = FOI_df,
aes(x = b_a, y = Mean_Presymp_Prop)) +
geom_point() +
rahul_man_figure_theme + theme_white_background
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Mean_Presymp_prop_of_FOI.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = FOI_df,
aes(x = b_a, y = Mean_Asymp_Prop)) +
geom_point() +
rahul_man_figure_theme + theme_white_background
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Mean_Asymp_prop_of_FOI.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = FOI_df,
aes(x = b_a, y = (R_0_S_1 + R_0_S_2)/R_0_NGM)) +
geom_point() + ylab("Prop Symp R_0") +
rahul_man_figure_theme + theme_white_background
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Symp_prop_of_R_0.png")
print(p)
dev.off## function (which = dev.cur())
## {
## if (which == 1)
## stop("cannot shut down device 1 (the null device)")
## .External(C_devoff, as.integer(which))
## dev.cur()
## }
## <bytecode: 0x7f955d574698>
## <environment: namespace:grDevices>p = ggplot(data = FOI_df,
aes(x = b_a, y = (R_0_P)/R_0_NGM)) +
geom_point() + ylab("Prop Pre-Symp R_0") +
rahul_man_figure_theme + theme_white_background
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Pre_Symp_prop_of_R_0.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = FOI_df,
aes(x = b_a, y = (R_0_A)/R_0_NGM)) +
geom_point() + ylab("Prop Asymp R_0") +
rahul_man_figure_theme + theme_white_background
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/b_a_vs_Asymp_prop_of_R_0.png")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data = FOI_df,
aes(x = b_a*b_p,
y = Mean_Symp_prop)) +
geom_point() +
rahul_man_figure_theme + theme_white_background
p
Figure 4 Alternate
Sub-funcitons
get_median_and_quantiles = function(single_b_a_value_FOI_data, type){
if(type == "Symptomatic"){
single_b_a_value_FOI_data = single_b_a_value_FOI_data %>%
mutate(FOI_contrib = lambda_FOI_S/lambda_FOI_total)
}else{
if(type == "Asymptomatic"){
single_b_a_value_FOI_data = single_b_a_value_FOI_data %>%
mutate(FOI_contrib = lambda_FOI_A/lambda_FOI_total)
}else{
if(type == "Presymptomatic"){
single_b_a_value_FOI_data = single_b_a_value_FOI_data %>%
mutate(FOI_contrib = lambda_FOI_P/lambda_FOI_total)
}else{
error("Invalid infection class")
}
}
}
single_b_a_value_summary = single_b_a_value_FOI_data %>%
group_by(time) %>%
summarize(median_FOI_contrib =
median(FOI_contrib),
upper_quantile_FOI_contrib = quantile(FOI_contrib, probs = 0.975),
lower_quantile_FOI_contrib = quantile(FOI_contrib, probs = 0.025)
) %>%
as.data.frame()%>%
mutate(Infection_Type = type)
}Define sampling times for FOI
head(Observed_data)## Y times obs_prop_positive
## 1 0 1 NA
## 2 0 2 0.00000000
## 3 2 3 0.25000000
## 4 2 4 0.05555556
## 5 7 5 0.15555556
## 6 0 6 0.00000000peak_time = Observed_data %>%
filter(Y == max(Y)) %>%
dplyr::select(times) %>%
as.numeric()
pre_peak_time = peak_time - 28
post_peak_time = peak_time + 28
p = ggplot(data = Observed_data, aes(x = times, y = Y)) + geom_point() + geom_line() + rahul_man_figure_theme +
geom_vline(xintercept = pre_peak_time, color = 'blue',
size = 1.5) +
geom_vline(xintercept = peak_time, color = 'red',
size = 1.5) +
geom_vline(xintercept = post_peak_time, color = 'orange',
size = 1.5) +
theme_white_background
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/sampling_points_for_FOI.png")
print(p)
dev.off()## quartz_off_screen
## 2Calculate summary statistics for FOI at each sampling time for each \(b_a\) including outlier value
b_a_points = unique(antibody_top_2_LL_from_b_a_profile_top_2_LL$b_a)
FOI_summary = data.frame(matrix(nrow = 0, ncol =6 ))
colnames(FOI_summary) = c("time" ,"median_FOI_contrib",
"upper_quantile_FOI_contrib",
"lower_quantile_FOI_contrib",
"Infection_Type", "b_a")
for(b_a_index in seq(1:length(b_a_points))){
print(b_a_index)
single_b_a_value = b_a_points[b_a_index]
param_combs_with_b_a_index_value =
antibody_top_2_LL_from_b_a_profile_top_2_LL %>%
filter(b_a == single_b_a_value)
pre_peak_single_b_a_value_FOI_data =
data.frame(matrix(nrow = 0, ncol = 13))
colnames(pre_peak_single_b_a_value_FOI_data) =
c("time",".id", "lambda_FOI_S" ,"lambda_FOI_A",
"lambda_FOI_P","lambda_FOI_total", "param_index")
post_peak_single_b_a_value_FOI_data =
data.frame(matrix(nrow = 0, ncol = 13))
colnames(post_peak_single_b_a_value_FOI_data) =
c("time",".id", "lambda_FOI_S" ,"lambda_FOI_A",
"lambda_FOI_P","lambda_FOI_total", "param_index")
peak_single_b_a_value_FOI_data =
data.frame(matrix(nrow = 0, ncol = 13))
colnames(peak_single_b_a_value_FOI_data) =
c("time",".id", "lambda_FOI_S" ,"lambda_FOI_A",
"lambda_FOI_P","lambda_FOI_total", "param_index")
for(param_index in
seq(1:nrow(param_combs_with_b_a_index_value))){
single_param_combo =
param_combs_with_b_a_index_value[param_index,]
param_combo_for_sim = single_param_combo %>%
dplyr::select(-msg, -iter_num, -param_index,
-loglik, -nfail, -trace_num, -loglist.se,
-Antibody_Mean_LL, -Antibody_LL_SE,
-Median_Herd_Immunity, -combo_num,
-sim_subset_index, -duration_of_symp_1,
-duration_of_symp_2, -duration_of_symp,
-gamma_total, -Beta,
-R_0_P, -R_0_A, -R_0_S_1, -R_0_S_2,
-R_0_NGM)
sim_data_single_param = simulate(nsim = 101,
seed = 12345,
times = Observed_data$times,
t0 = t0,
rprocess = pomp::euler(rproc,delta.t = 1),
params = param_combo_for_sim,
paramnames = paramnames,
statenames = statenames,
obsnames = obsnames,
accumvars = acumvarnames,
rinit = init,
rmeas = rmeas,
partrans = par_trans,
covar = covar,
format = "data.frame")
#head(sim_data)
sim_data_single_param_FOI_data = sim_data_single_param %>%
dplyr::select(time, .id, I_P, I_S_1, I_S_2,
I_A, N, beta_t)
sim_data_single_param_FOI_data = sim_data_single_param_FOI_data %>%
mutate(lambda_FOI_S = (beta_t*(I_S_1 + I_S_2))/N,
lambda_FOI_A = (beta_t*param_combo_for_sim$b_a*I_A)/N,
lambda_FOI_P = (beta_t*param_combo_for_sim$b_p*I_P)/N) %>%
mutate(lambda_FOI_total =
lambda_FOI_S + lambda_FOI_A + lambda_FOI_P) %>%
dplyr::select(time, .id, lambda_FOI_S, lambda_FOI_A, lambda_FOI_P,
lambda_FOI_total) %>%
mutate(param_index = param_index)
pre_peak_single_param_FOI_data =
sim_data_single_param_FOI_data %>%
filter(time == pre_peak_time)
peak_time_single_param_FOI_data =
sim_data_single_param_FOI_data %>%
filter(time == peak_time)
post_peak_single_param_FOI_data =
sim_data_single_param_FOI_data %>%
filter(time == post_peak_time)
pre_peak_single_b_a_value_FOI_data = rbind(
pre_peak_single_b_a_value_FOI_data,
pre_peak_single_param_FOI_data
)
peak_single_b_a_value_FOI_data = rbind(
peak_single_b_a_value_FOI_data,
peak_time_single_param_FOI_data
)
post_peak_single_b_a_value_FOI_data = rbind(
post_peak_single_b_a_value_FOI_data,
post_peak_single_param_FOI_data
)
}
### Calculate % contribution of FOI from each category
single_b_a_value_FOI_data = rbind(
pre_peak_single_b_a_value_FOI_data,
peak_single_b_a_value_FOI_data) %>%
rbind(post_peak_single_b_a_value_FOI_data)
single_b_a_Presymp_summary = get_median_and_quantiles(
single_b_a_value_FOI_data = single_b_a_value_FOI_data,
type = "Presymptomatic")
single_b_a_Symp_summary = get_median_and_quantiles(
single_b_a_value_FOI_data = single_b_a_value_FOI_data,
type = "Symptomatic")
single_b_a_Asymp_summary = get_median_and_quantiles(
single_b_a_value_FOI_data = single_b_a_value_FOI_data,
type = "Asymptomatic")
single_b_a_FOI_summary = rbind(single_b_a_Presymp_summary,
single_b_a_Symp_summary) %>%
rbind(single_b_a_Asymp_summary) %>%
mutate(b_a = single_b_a_value)
FOI_summary = FOI_summary %>%
rbind(single_b_a_FOI_summary)
}## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21Isolate times
peak_FOI_summary = FOI_summary %>%
filter(time == peak_time)
pre_peak_FOI_summary = FOI_summary %>%
filter(time == pre_peak_time)
post_peak_FOI_summary = FOI_summary %>%
filter(time == post_peak_time)Plot FOI
Peak
p = ggplot(data= peak_FOI_summary,
aes(x = b_a,
y = median_FOI_contrib,
color = Infection_Type)) +
geom_point() +
geom_line() + geom_ribbon(aes(ymin = lower_quantile_FOI_contrib,
ymax = upper_quantile_FOI_contrib,
fill = Infection_Type), alpha = 0.5) +
rahul_man_figure_theme +
theme_white_background + xlab("Asymptomatic Transmission Rate ($b_a$) ") +
ylab("% Contribution to the Force of Infection at Peak")
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/b_a_vs_FOI_contrib_line_plot_peak_include_outlier.png")
print(p)
dev.off()## quartz_off_screen
## 2Pre-Peak
p = ggplot(data= pre_peak_FOI_summary,
aes(x = b_a,
y = median_FOI_contrib,
color = Infection_Type)) +
geom_point() +
geom_line() + geom_ribbon(aes(ymin = lower_quantile_FOI_contrib,
ymax = upper_quantile_FOI_contrib,
fill = Infection_Type), alpha = 0.5) +
rahul_man_figure_theme +
theme_white_background + xlab("Asymptomatic Transmission Rate ($b_a$) ") +
ylab("% Contribution to the Force of Infection 4 Weeks Before Peak")
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/b_a_vs_FOI_contrib_line_plot_pre_peak_include_outlier.png")
print(p)
dev.off()## quartz_off_screen
## 2Post-Peak
p = ggplot(data= post_peak_FOI_summary,
aes(x = b_a,
y = median_FOI_contrib,
color = Infection_Type)) +
geom_point() +
geom_line() + geom_ribbon(aes(ymin = lower_quantile_FOI_contrib,
ymax = upper_quantile_FOI_contrib,
fill = Infection_Type), alpha = 0.5) +
rahul_man_figure_theme +
theme_white_background + xlab("Asymptomatic Transmission Rate ($b_a$) ") +
ylab("% Contribution to the Force of Infection 4 Weeks After Peak")
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/b_a_vs_FOI_contrib_line_plot_post_peak_include_outlier.png")
print(p)
dev.off()## quartz_off_screen
## 2Calculate summary statistics for FOI at each sampling time for each \(b_a\) excluding outlier value
second_outlier = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_outlier %>%
filter(b_p < 0.02)
b_a_points_no_outlier = unique(antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers$b_a)
FOI_summary_no_outlier = data.frame(matrix(nrow = 0, ncol =6 ))
colnames(FOI_summary_no_outlier) = c("time" ,"median_FOI_contrib",
"upper_quantile_FOI_contrib",
"lower_quantile_FOI_contrib",
"Infection_Type", "b_a")
for(b_a_index in seq(1:length(b_a_points_no_outlier))){
print(b_a_index)
single_b_a_value = b_a_points_no_outlier[b_a_index]
param_combs_with_b_a_index_value_no_outlier =
antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers %>%
filter(b_a == single_b_a_value)
pre_peak_single_b_a_value_FOI_data =
data.frame(matrix(nrow = 0, ncol = 13))
colnames(pre_peak_single_b_a_value_FOI_data) =
c("time",".id", "lambda_FOI_S" ,"lambda_FOI_A",
"lambda_FOI_P","lambda_FOI_total", "param_index")
post_peak_single_b_a_value_FOI_data =
data.frame(matrix(nrow = 0, ncol = 13))
colnames(post_peak_single_b_a_value_FOI_data) =
c("time",".id", "lambda_FOI_S" ,"lambda_FOI_A",
"lambda_FOI_P","lambda_FOI_total", "param_index")
peak_single_b_a_value_FOI_data =
data.frame(matrix(nrow = 0, ncol = 13))
colnames(peak_single_b_a_value_FOI_data) =
c("time",".id", "lambda_FOI_S" ,"lambda_FOI_A",
"lambda_FOI_P","lambda_FOI_total", "param_index")
for(param_index in
seq(1:nrow(param_combs_with_b_a_index_value_no_outlier))){
single_param_combo =
param_combs_with_b_a_index_value_no_outlier[param_index,]
param_combo_for_sim = single_param_combo %>%
dplyr::select(-msg, -iter_num, -param_index,
-loglik, -nfail, -trace_num, -loglist.se,
-Antibody_Mean_LL, -Antibody_LL_SE,
-Median_Herd_Immunity, -combo_num,
-sim_subset_index, -duration_of_symp_1,
-duration_of_symp_2, -duration_of_symp,
-gamma_total, -Beta,
-R_0_P, -R_0_A, -R_0_S_1, -R_0_S_2,
-R_0_NGM)
sim_data_single_param = simulate(nsim = 101,
seed = 12345,
times = Observed_data$times,
t0 = t0,
rprocess = pomp::euler(rproc,delta.t = 1),
params = param_combo_for_sim,
paramnames = paramnames,
statenames = statenames,
obsnames = obsnames,
accumvars = acumvarnames,
rinit = init,
rmeas = rmeas,
partrans = par_trans,
covar = covar,
format = "data.frame")
#head(sim_data)
sim_data_single_param_FOI_data = sim_data_single_param %>%
dplyr::select(time, .id, I_P, I_S_1, I_S_2,
I_A, N, beta_t)
sim_data_single_param_FOI_data = sim_data_single_param_FOI_data %>%
mutate(lambda_FOI_S = (beta_t*(I_S_1 + I_S_2))/N,
lambda_FOI_A = (beta_t*param_combo_for_sim$b_a*I_A)/N,
lambda_FOI_P = (beta_t*param_combo_for_sim$b_p*I_P)/N) %>%
mutate(lambda_FOI_total =
lambda_FOI_S + lambda_FOI_A + lambda_FOI_P) %>%
dplyr::select(time, .id, lambda_FOI_S, lambda_FOI_A, lambda_FOI_P,
lambda_FOI_total) %>%
mutate(param_index = param_index)
pre_peak_single_param_FOI_data =
sim_data_single_param_FOI_data %>%
filter(time == pre_peak_time)
peak_time_single_param_FOI_data =
sim_data_single_param_FOI_data %>%
filter(time == peak_time)
post_peak_single_param_FOI_data =
sim_data_single_param_FOI_data %>%
filter(time == post_peak_time)
pre_peak_single_b_a_value_FOI_data = rbind(
pre_peak_single_b_a_value_FOI_data,
pre_peak_single_param_FOI_data
)
peak_single_b_a_value_FOI_data = rbind(
peak_single_b_a_value_FOI_data,
peak_time_single_param_FOI_data
)
post_peak_single_b_a_value_FOI_data = rbind(
post_peak_single_b_a_value_FOI_data,
post_peak_single_param_FOI_data
)
}
### Calculate % contribution of FOI from each category
single_b_a_value_FOI_data = rbind(
pre_peak_single_b_a_value_FOI_data,
peak_single_b_a_value_FOI_data) %>%
rbind(post_peak_single_b_a_value_FOI_data)
single_b_a_Presymp_summary = get_median_and_quantiles(
single_b_a_value_FOI_data = single_b_a_value_FOI_data,
type = "Presymptomatic")
single_b_a_Symp_summary = get_median_and_quantiles(
single_b_a_value_FOI_data = single_b_a_value_FOI_data,
type = "Symptomatic")
single_b_a_Asymp_summary = get_median_and_quantiles(
single_b_a_value_FOI_data = single_b_a_value_FOI_data,
type = "Asymptomatic")
single_b_a_FOI_summary = rbind(single_b_a_Presymp_summary,
single_b_a_Symp_summary) %>%
rbind(single_b_a_Asymp_summary) %>%
mutate(b_a = single_b_a_value)
FOI_summary_no_outlier = FOI_summary_no_outlier %>%
rbind(single_b_a_FOI_summary)
}## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21Isolate times
peak_FOI_summary_no_outlier = FOI_summary_no_outlier %>%
filter(time == peak_time)
pre_peak_FOI_summary_no_outlier = FOI_summary_no_outlier %>%
filter(time == pre_peak_time)
post_peak_FOI_summary_no_outlier = FOI_summary_no_outlier %>%
filter(time == post_peak_time)Peak
p = ggplot(data= peak_FOI_summary_no_outlier,
aes(x = b_a,
y = median_FOI_contrib,
color = Infection_Type)) +
geom_point() +
geom_line() + geom_ribbon(aes(ymin = lower_quantile_FOI_contrib,
ymax = upper_quantile_FOI_contrib,
fill = Infection_Type), alpha = 0.5) +
rahul_man_figure_theme +
theme_white_background + xlab("Asymptomatic Transmission Rate ($b_a$) ") +
ylab("% Contribution to the Force of Infection at Peak")
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/b_a_vs_FOI_contrib_line_plot_peak_no_outlier.png")
print(p)
dev.off()## quartz_off_screen
## 2Pre-Peak
p = ggplot(data= pre_peak_FOI_summary_no_outlier,
aes(x = b_a,
y = median_FOI_contrib,
color = Infection_Type)) +
geom_point() +
geom_line() + geom_ribbon(aes(ymin = lower_quantile_FOI_contrib,
ymax = upper_quantile_FOI_contrib,
fill = Infection_Type), alpha = 0.5) +
rahul_man_figure_theme +
theme_white_background + xlab("Asymptomatic Transmission Rate ($b_a$) ") +
ylab("% Contribution to the Force of Infection 4 Weeks Before Peak")
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/b_a_vs_FOI_contrib_line_plot_pre_peak_no_outlier.png")
print(p)
dev.off()## quartz_off_screen
## 2Post-Peak
p = ggplot(data= post_peak_FOI_summary_no_outlier,
aes(x = b_a,
y = median_FOI_contrib,
color = Infection_Type)) +
geom_point() +
geom_line() + geom_ribbon(aes(ymin = lower_quantile_FOI_contrib,
ymax = upper_quantile_FOI_contrib,
fill = Infection_Type), alpha = 0.5) +
rahul_man_figure_theme +
theme_white_background + xlab("Asymptomatic Transmission Rate ($b_a$) ") +
ylab("% Contribution to the Force of Infection 4 Weeks After Peak")
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/b_a_vs_FOI_contrib_line_plot_post_peak_no_outlier.png")
print(p)
dev.off()## quartz_off_screen
## 2Peak median density plot
p = ggplot(data=peak_FOI_summary_no_outlier,
aes(x=factor(b_a), y=median_FOI_contrib, fill=Infection_Type)) +
geom_bar(stat="identity", position="stack") +
geom_errorbar(aes(ymin = lower_quantile_FOI_contrib,
ymax = upper_quantile_FOI_contrib))
p
### Sub-function for stacking data
library(tidyr)## Warning: package 'tidyr' was built under R version 3.5.2##
## Attaching package: 'tidyr'## The following object is masked from 'package:reshape2':
##
## smithsstack_FOI_error_bars = function(FOI_summary_no_outlier){
test_df_asymp =FOI_summary_no_outlier %>%
filter(Infection_Type == "Asymptomatic")
spreaded_median_FOI = FOI_summary_no_outlier %>%
spread(Infection_Type, median_FOI_contrib)
median_asymp_FOI = spreaded_median_FOI %>%
dplyr::select(b_a, Asymptomatic_Median = Asymptomatic) %>%
na.omit()
median_symp_FOI = spreaded_median_FOI %>%
dplyr::select(b_a, Symptomatic_Median = Symptomatic) %>%
na.omit()
median_presymp_FOI = spreaded_median_FOI %>%
dplyr::select(b_a, Presymptomatic_Median = Presymptomatic) %>%
na.omit()
median_combined_FOI = median_asymp_FOI %>%
join(median_presymp_FOI) %>%
join(median_symp_FOI)
FOI_summary_no_outlier_with_eb = FOI_summary_no_outlier %>%
join(median_combined_FOI) %>%
mutate(Add_Presymp = as.numeric(Infection_Type == "Asymptomatic"),
Add_Symp = as.numeric(Infection_Type %in% c("Asymptomatic", "Presymptomatic" )),
Stacked_Median_Adjustment = Add_Presymp*Presymptomatic_Median +
Add_Symp*Symptomatic_Median)
FOI_summary_no_outlier_with_eb = FOI_summary_no_outlier_with_eb%>%
mutate(b_a_factor = as.factor(b_a))
return(FOI_summary_no_outlier_with_eb)
}Stack error bars
peak_FOI_summary_no_outlier_with_eb = stack_FOI_error_bars(
FOI_summary_no_outlier = peak_FOI_summary_no_outlier)## Joining by: b_a
## Joining by: b_a
## Joining by: b_apre_peak_FOI_summary_no_outlier_with_eb = stack_FOI_error_bars(
FOI_summary_no_outlier = pre_peak_FOI_summary_no_outlier)## Joining by: b_a
## Joining by: b_a
## Joining by: b_apost_peak_FOI_summary_no_outlier_with_eb = stack_FOI_error_bars(
FOI_summary_no_outlier = post_peak_FOI_summary_no_outlier
)## Joining by: b_a
## Joining by: b_a
## Joining by: b_aall_b_a_levels = levels(peak_FOI_summary_no_outlier_with_eb$b_a_factor)
b_a_level_tickmarks = all_b_a_levels[c(1,7,13,18,21)]
b_a_tick_labels = c("0", "0.24", "0.52", "0.79", "0.97")Plot stacked bar plot for peak
p = ggplot(data=peak_FOI_summary_no_outlier_with_eb,
aes(x=b_a_factor, y=median_FOI_contrib, fill=Infection_Type)) +
geom_bar(stat="identity", position="stack") +
geom_errorbar(aes(ymin = Stacked_Median_Adjustment + lower_quantile_FOI_contrib,
ymax = Stacked_Median_Adjustment + upper_quantile_FOI_contrib)) +
scale_x_discrete(breaks = b_a_level_tickmarks,
labels = b_a_tick_labels) +
xlab('Relative Asymptomatic \n Transmission Rate') +
ylab("Contribution to Force of Infection At Peak ") +
scale_fill_discrete(name = "Infection \n Type") +
rahul_man_figure_theme + theme_white_background +
theme(axis.title.x = element_text(face = "plain", size = 12),
axis.title.y = element_text(face = "plain", size = 12),
axis.text.x = element_text(face = "plain", size = 11),
axis.text.y = element_text(face = "plain", size = 11),
legend.text = element_text(face = "plain", size = 11),
legend.title = element_text(face = "plain", size = 12))+
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/FOI_Stacked_Bar_Plot_Peak_exclude_two_outliers.png")
print(p)
dev.off()## quartz_off_screen
## 2pdf("../Figures/Profiles/N_12_Model/Man_Figs/stacked_bar_plot_b_a_vs_cont_to_FOI_at_peak_exclude_two_outliers.pdf")
print(p)
dev.off()## quartz_off_screen
## 2p = ggplot(data=peak_FOI_summary_no_outlier_with_eb,
aes(x=b_a_factor, y=median_FOI_contrib, fill=Infection_Type)) +
geom_bar(stat="identity", position="stack") +
geom_errorbar(aes(ymin = Stacked_Median_Adjustment + lower_quantile_FOI_contrib,
ymax = Stacked_Median_Adjustment + upper_quantile_FOI_contrib)) +
scale_x_discrete(breaks = b_a_level_tickmarks,
labels = b_a_tick_labels) +
xlab('Relative Asymptomatic \n Transmission Rate') +
ylab("Contribution to Force of Infection At Peak ") +
scale_fill_discrete(name = "Infection \n Type") +
rahul_man_figure_theme + theme_white_background +
theme(axis.title.x = element_text(face = "plain", size = 12),
axis.title.y = element_text(face = "plain", size = 12),
axis.text.x = element_text(face = "plain", size = 11),
axis.text.y = element_text(face = "plain", size = 11),
legend.text = element_text(face = "plain", size = 11),
legend.title = element_text(face = "plain", size = 12)) +
scale_fill_grey(name = "Infection \n Type")+
theme(axis.line = element_line(colour = 'black', size = 1))+
theme(axis.ticks = element_line(colour = "black", size = 1.5)) ## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.p
pdf("../Figures/Profiles/N_12_Model/Man_Figs/stacked_bar_plot_b_a_vs_cont_to_FOI_at_peak_exclude_two_outliers_greyscale.pdf")
print(p)
dev.off()## quartz_off_screen
## 2Plot stacked bar plot for 4 weeks before peak
p = ggplot(data=pre_peak_FOI_summary_no_outlier_with_eb,
aes(x=b_a_factor, y=median_FOI_contrib, fill=Infection_Type)) +
geom_bar(stat="identity", position="stack") +
geom_errorbar(aes(ymin = Stacked_Median_Adjustment + lower_quantile_FOI_contrib,
ymax = Stacked_Median_Adjustment + upper_quantile_FOI_contrib)) +
scale_x_discrete(breaks = b_a_level_tickmarks,
labels = b_a_tick_labels) +
xlab('Relative Asymptomatic \n Transmission Rate') +
ylab("Contribution to Force of Infection \n 4 Weeks Before Peak ") +
scale_fill_discrete(name = "Infection \n Type") +
rahul_man_figure_theme + theme_white_background +
theme(axis.title.x = element_text(face = "plain", size = 12),
axis.title.y = element_text(face = "plain", size = 12),
axis.text.x = element_text(face = "plain", size = 11),
axis.text.y = element_text(face = "plain", size = 11),
legend.text = element_text(face = "plain", size = 11),
legend.title = element_text(face = "plain", size = 12))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/FOI_Stacked_Bar_Plot_Pre_Peak_exclude_two_outliers.png")
print(p)
dev.off()## quartz_off_screen
## 2Plot stacked bar plot for 4 weeks after peak
p = ggplot(data=post_peak_FOI_summary_no_outlier_with_eb,
aes(x=b_a_factor, y=median_FOI_contrib, fill=Infection_Type)) +
geom_bar(stat="identity", position="stack") +
geom_errorbar(aes(ymin = Stacked_Median_Adjustment + lower_quantile_FOI_contrib,
ymax = Stacked_Median_Adjustment + upper_quantile_FOI_contrib)) +
scale_x_discrete(breaks = b_a_level_tickmarks,
labels = b_a_tick_labels) +
xlab('Relative Asymptomatic \n Transmission Rate') +
ylab("Contribution to Force of Infection \n 4 Weeks After Peak ") +
scale_fill_discrete(name = "Infection \n Type") +
rahul_man_figure_theme + theme_white_background +
theme(axis.title.x = element_text(face = "plain", size = 12),
axis.title.y = element_text(face = "plain", size = 12),
axis.text.x = element_text(face = "plain", size = 11),
axis.text.y = element_text(face = "plain", size = 11),
legend.text = element_text(face = "plain", size = 11),
legend.title = element_text(face = "plain", size = 12))
p
png("../Figures/Profiles/N_12_Model/top_2_LL_sim_plots/FOI_plots/Timepoint_Analysis/FOI_Stacked_Bar_Plot_Post_Peak_exclude_two_outliers.png")
print(p)
dev.off()## quartz_off_screen
## 2Sup Figures
p = ggplot(data = antibody_top_2_LL_from_b_a_profile_top_2_LL_no_2_outliers,
aes(x = R_0,
y = R_0_NGM,
color = b_p)) +
geom_point(size = 5) +
scale_color_viridis_c() +
rahul_man_figure_theme +
theme_white_background +
scale_x_continuous(breaks=c(seq(2,10,1), 15, 18)) +
scale_y_continuous(breaks=seq(2,5,1)) +
coord_cartesian(expand = FALSE, #turn off axis expansion (padding)
xlim = c(1.75, 9), ylim = c(1.75, 5.25)) #manually set limits
p
png(paste0("../Figures/Profiles/", model_name, "_Model/Sup_Figs/",
"R_0_vs_R_0_NGM_color_by_b_p", model_name,
"_model_top_2_antibody_LL_from_b_a_profile_peak_LL_no_2_outliers.png"))
print(p)
dev.off()## quartz_off_screen
## 2